Inversion of Complex Body Waves

A numerical method to deconvolve complex body waves into a multiple shock sequence is developed. With the assumption that all the constituent events of a multiple shock have identical fault geometry and depth, the far-field source time function is obtained as a superposiUon of ramp functions. The height and the onset time of the ramp functions are determined by matching the synthetic waveforms with the observed ones in the least-square sense. The individual events are then identified by pairs of ramp functions or discrete trapezoidal pulses in the source time sequence. The method can be used for the analysis of both single and multi-station data. Teleseismic long-period P waves from the 1976 Guatemala earthquake are analyzed as a test of our method. The present method provides a useful tool for a systematic analysis of multiple event sequences. INTRODUCTION The spectra and waveforms of seismic body waves provide important information on the details of the source rupture process. In frequency domain analysis, the lowfrequency asymptote and the corner frequency of the displacement spectrum are used to estimate the seismic moment and the source dimension (Brune, 1970). In time-domain analysis of body waves, the observed waveforms are modeled by a source time function, and the time constants associated with it are interpreted in terms of the source dimension and the particle velocity of the fault motion (Aki, 1968; Haskell, 1969; Kanamori, 1972; Abe, 1974). When the observed body waveforms are relatively simple, the modeling can be done by using either forward or inverse methods. Langston (1976) and Burdick and Mellman {1976) used a time-domain inversion method to determine some of the complexities of the source time function. For a very large earthquake, however, the body waveform is extremely complex, and no standard method is available for the inversion. Several attempts have been made to explain the complexity of body waves from large earthquakes by using a multiple event model. Earlier attempts consisted of identifying distinct phases in the body-wave signal as the beginning of distinct events and locating them with respect to the first one (Imamura, 1937, p. 267; Miyamura et al., 1964; Wyss and Brune, 1967). In more recent studies, synthetic seismograms were used to make a more quantitative interpretation (Fukao, 1972; Chung and Kanamori, 1976). Kanamori and Stewart (1978) modeled the waveforms of P waves from the 1976 Guatemala earthquake by matching them, in the leastsquare sense, with a series of trapezoidal source time functions. Rial (1978) modeled the Caracas earthquake by using three distinct events. Boatwright (1980) employed a direct inversion of the body waves from the 1979 St. Elias, Alaska, earthquake to resolve a few subevents. The complexity of the source time function reflects the heterogeneity in the mechanical properties in the fault zone, which is often characterized by asperities or barriers. Many recent studies have suggested the importance of asperities in various seismological problems, such as the nature of strong ground motion (Das and Aki, * Present address. Physms Department, Yokohama City University, Yokohama 236, Japan. 491 4 9 2 MASAYUKI KIKUCHI AND HIROO KANAMORI 1977; Ebel and Helmberger, 1981), foreshocks (Jones and Molnar, 1979), seismicity patterns (Wesson and Ellsworth, 1973; Kanamori, 1981), and the regional variation of rupture mode (Lay and Kanamori, 1981). In view of this importance, it is desirable to develop a systematic method for inversion of complex body waves consisting of contributions from several sources. This is obviously a difficult problem. For example, the 1976 Guatemala earthquake was modeled by about 10 pulses, each representing a distinct seismic source. Even in the simplest source model, about six parameters [the seismic moment, three fault parameters, two time constants (e.g., rise time and pulse width)] are necessary to describe each source. Thus, if we are to model a multiple shock with 10 distinct events, about 100 parameters, including the origin time and the coordinates of the individual events, would have to be determined. In view of the amount, the quality and the limited bandwidth of the data usually available for this type of modeling, it would be very difficult to determine all of the parameters. Furthermore, in view of the complexity of the structure near the source, along the path, and near the receiver, it would not be easy to prove that all the complexities in the body wave form are due to the source. Because of these difficulties, we will be mainly concerned with the gross complexities of multiple events rather than with the minute details of the source function, and a number of simplifications will be made. Inevitably, a certain amount of nonuniqueness and subjectivity exists. The validity of the model should eventually be judged in the light of other data such as local strong-motion data, distribution and geometry of surface breaks, and macroseismic data. As we will show in the later sections, the method we present here can invert complex observed seismograms into a source time function in a systematic and reasonably objective way, thereby providing a means for interpreting complex observed records in terms of asperities and barriers in the fault zone. M E T H O D In an infinite homogeneous space, the far-field body waveform due to a shear dislocation source is given by [e.g., (10) in Haskell, 1964] Re# f f D($, t r/c) dA Uc(f, t) 4~rpc3ro (1) where A = dislocation surface, ~ = a variable point on A, 2 = an observation point, t = time, r = ]2 ~1, r0 = the average of r,/9(~, t) = relative slip velocity, Re = radiation pattern, p = density,/~ = rigidity, and c = body-wave velocity. When the source region is small, the travel time r/c in (1) can be approximated by its average, ro/c. The waveform is then given by Re Ue(2, t) 4wpc3r------° S(t ro/c) (2) where S(t) is the far-field source time function defined by s ( t ) = ~ f fAD(g , t )dA. (3) Here we assume that the time history of dislocation at a point is given by a INVERSION OF COMPLEX BODY WAVES 493 function of the time measured from the arrival of a rupture front. Let t'(~) be the arrival time at a point ~, then the dislocation function is expressed as D(~ , t) = D ( t t ' (~ ) ) . (4) Noting that d A = ( d A / d t ' ) d t ' is the area swept by the rupture front during the time interval d t ' , we can write equation (3) as S ( t ) = it D ( t t ' ) A ( t ' ) d t ' (5) where dot denotes the time derivative. Thus the far-field source time function is expressed by a convolution of the dislocation velocity and the fault area expansion rate. We assume that the disolcation time history is given by a ramp function with rise time r as D ( t ) = D o s ( t ) where Do is the final dislocation and s (t) is the unit ramp function s ( t ) = • 0 ¢ characterized by abrupt changes of the fault area If rupture propagation is expansion rate, then A ( t ) = ~ A A , H ( t t~) (6) l where ha t is the increment of the fault area expansion rate at time t,, a n d H ( t ) is the Heaviside step function. The source time function S ( t ) is then given by superposition of ramp functions where S ( t ) = ~ m , s ( t t ,) (7) rn~ = It D o A A , . For example, in case of unilateral rupture propagation, A = W v t ( O <~ t <~ T ) (W = fault width, v -rupture velocity) and the far-field source time function can be described by a pair of pulses as follows m l = It W v D o , tl --O, m2 = --it W v D o , t2 = T. 494 MASAYUKI KIKUCHI AND I-IIROO KANAMORI When these pulses are convolved with s( t ) , a trapezoidal far-field source time function is produced. In case of asymmetric bilateral rupture, A = 2 w v t for O <~ t < T~ w v t for T~ < t <<T2